Fluctuating surface-current formulation of radiative heat transfer for arbitrary geometries 



Alejandro W. Rodriguez, 1 ' 2 M. T. H. Reid, 2 and Steven G. Johnson 2 

School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138 
2 Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139 

We describe a fluctuating surface-current formulation of radiative heat transfer, applicable to arbitrary geome- 
tries, that directly exploits standard, efficient, and sophisticated techniques from the boundary-element method. 
We validate as well as extend previous results for spheres and cylinders, and also compute the heat transfer in 
a more complicated geometry consisting of two interlocked rings. Finally, we demonstrate that the method can 
be readily adapted to compute the spatial distribution of heat flux on the surface of the interacting bodies. 

PACS numbers: 



Quantum and thermal fluctuations of charges in otherwise 
neutral bodies lead to stochastic electromagnetic (EM) fields 
everywhere in space. In non-equilibrium situations involv- 
ing bodies at different temperatures, these fields mediate en- 
ergy exchange from the hotter to the colder bodies, a process 
known as radiative heat transfer. Although the basic theoret- 
ical formalism for studying heat transfer was laid out decades 
ago [1-6], only recently have experiments reached the preci- 
sion required to measure them at the microscale [7, 8], spark- 
ing renewed interest in the study of these interactions in com- 
plex geometries that deviate from the simple parallel-plate 
structures of the past. In this letter, we propose a novel formu- 
lation of radiative heat transfer for arbitrary geometries that is 
based on the fluctuating surface-current (FSC) method of clas- 
sical EM fields [9]. Unlike previous scattering formulations 
based on basis expansions of the field unknowns best suited to 
special [10-14] or non-interleaved periodic [15] geometries, 
or formulations based on expensive, brute-force time-domain 
simulations [16], this approach allows direct application of the 
boundary element method (BEM): a mature and sophisticated 
surface-integral equation (SIE) formulation of the scattering 
problem in which the EM fields are determined by the solu- 
tion of an algebraic equation involving a smaller set of surface 
unknowns (fictitious surface currents in the surfaces of the ob- 
jects [17]). In what follows, we briefly review the SIE method, 
derive an FSC equation for the heat transfer between two bod- 
ies, and demonstrate its correctness by checking it against (as 
well as extending) previous results for spheres and cylinders. 
To demonstrate the generality of this method, we compute the 
heat transfer in a complicated geometry that lies beyond the 
reach of other formulations, as well as show that it can be 
readily adapted to obtain the spatial distribution of flux pat- 
tern at the surface of the bodies. 

The radiative heat transfer between two objects 1 and 2 at 
local temperatures T 1 and T 2 can be written as [5, 6]: 
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where 9(cj,T) = Huj/[exp(huj/kBT) — 1] is the Planck en- 
ergy per oscillator at temperature T, and is an ensemble- 
averaged flux spectrum into object 2 due to random currents 
in object 1 (defined more precisely below via the fluctuation- 
dissipation theorem [1, 18]). The only question is how to com- 



pute <£, which naively involves a cumbersome number of scat- 
tering calculations. 

Formulation: We begin by presenting our final result for 
<£, which is derived and validated below. Consider homoge- 
neous objects 1 and 2 separated by a lossless medium 0. Let 
T r denote the 6 x 6 Green's function T r (x, y) = T r (x - y) 
of the homogeneous medium r at a given uj (known analyt- 
ically [19]), relating 6-component electric (J) and magnetic 
(M) currents £ = (J;M) [";" denoting vertical concatena- 
tion] to 6-component electric (E) and magnetic (H) fields 
<Mx) = (E; H) = r r * £ = /d 3 yr(x,y)e(y) via a con- 
volution (*). Remarkably, we find that can be expressed 
purely in terms of interactions of fictitious surface currents 
located on the interfaces of the objects. Let be a basis 
of 6-component tangential vector fields on the surface of ob- 
ject r, so that any surface current £ r can be written in the form 
T(x) = E n <SW for coefficients x r n . In BEM, f3 n is 
typically a piecewise-polynomial "element" function defined 
within discretized patches of each surface [17]. However, one 
could just as easily choose f3 n to be a spherical harmonic or 
some other "spectral" Fourier-like basis [13]. The key point 
is that f3 n is an arbitrary basis of surface vector fields; unlike 
scattering-matrix formulations [11-13], it need not consist of 
"incoming" or "outgoing" waves nor satisfy any wave equa- 
tion. Our final result is the compact expression: 
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with symG = \{G + G*), where * denotes conjugate- 
transpose. The G and W matrices relate surface currents f3 n 
to surface-tangential fields T * /3 m or vice versa. Specifically, 
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where (ip,(/))r — <ffi> r ty*<\> * s the standard inner product over 
the surface of medium r (over both surfaces and both sets of 



2 



basis functions if r = 0), and 



W 11 
W 21 



W 12 
W 22 



w 



G° 



G 1 



+ 



G 2 



G 1 



G 2 



(4) 

is the BEM matrix inverse, used to solve SIE scattering 
problems as reviewed below, which relates incident fields to 
"equivalent" surface currents. In particular, W 21 relates inci- 
dent fields at the surface of object 2 to the equivalent currents 
at the surface of object 1 . Equation (2) is computationally con- 
venient because it only involves standard matrices that arise in 
BEM calculations [17], with no explicit need for evaluation of 
fields or sources in the volumes, separation of incoming and 
outgoing waves, integration of Poynting fluxes, or any addi- 
tional scattering calculations. As explained below, one can 
also obtain spatially resolved Poynting fluxes on the surfaces 
of the objects, as well as the emissivity of a single object, by 
a slight modification of Eq. (2). 

In addition to its computational elegance, Eq. (2) alge- 
braically captures crucial physical properties of <£. The stan- 
dard definiteness properties of the Green's functions (currents 
do nonnegative work) imply that sym G r is negative semidef- 
inite and hence it has a Cholesky factorization symG r = 
— U r *U r where U r is upper-triangular. It follows that $ = 
^Tr[Z*Z] = ±\\Z\\ 2 F where Z = U 2 W 21 U X \ is a 
weighted Frobenius norm of the interaction matrix W 21 , and 
hence > as required. Furthermore, reciprocity (sym- 
metry of <1> under 1 o 2 interchange) corresponds to sim- 
ple symmetries of the matrices. Inspection of T shows that 
r(y,x) T = Sr(x,y)S, where S = S T = S' 1 = <S* is the 
matrix that flips the sign of the magnetic components, and it 
follows from (3) that G T = SGS and W T = SWS where 
S = S T = 5 -1 = S* is the matrix that flips the signs of the 
magnetic basis coefficients and swaps the coefficients of f3 n 
and f3 n . (For convenience, we assume f3 n to be real, which is 
true in the case of RWG basis functions [17].) It follows that 
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where the S factors cancel, leading to the 1 ^ 2 exchange. 

Derivation: The key to our derivation of (2) is the SIE 
formulation of EM scattering [17, 20], which we briefly re- 
view here. Consider the fields (jf = (/> r+ + (j) r ~ in each 
region r, where r+ is the "incident" field due to sources 
within medium r, and (j) r ~ is the "scattered" field due to 
both interface reflections and sources in the other media. The 
core idea in the SIE formulation is the principle of equiva- 
lence [20], which states that the scattered field (jf~ can be 
expressed as the field of some fictitious electric and magnetic 
surface currents £ r located on the boundary of region r, act- 
ing within an infinite homogeneous medium r. In particular, 
the field in is = T° ★ (f 1 + £ 2 ). Remarkably, 



the same currents with a sign flip describe scattered fields 
in the interiors of the two objects [20]: (j) r ~ = — T r ★ £ r 
for r = 1,2. These currents £ r , in turn, are completely de- 
termined by the boundary condition of continuous tangen- 
tial fields = r | r at the r = 1,2 interfaces, giving 
theSIEs (r° + r r )*«f + r°*«f- r | r = r +-0°+| r for 
£ r in terms of the incident fields. To obtain a discrete set 
of equations, one expands £ r = x r n fi r n in a basis fi r n as 
above, and then takes the inner product of both sides of the 
SIEs with (a Galerkin discretization) to obtain a matrix 
"BEM" equation W~ 1 x = s in terms of exactly the W ma- 
trix from Eq. (4), current coefficients x = (x l ]x 2 ), and a 
right-hand "source" term s = (s 1 ; s 2 ) from the incident fields: 

^ = (^,0 r+ -0° + >r[17]. 

To compute we start by considering the flux <£ s into ob- 
ject 2 due to a single dipole source a 1 within object 1, so 
that 1+ = T 1 * a 1 and all other incident fields are zero. 
This corresponds to a right-hand side s = (s^O) where 
s m = (At^r 1 ★ a 1 )! in the BEM equation. <£ s is the re- 
sulting absorbed power in object 2, equal to the net incom- 
ing Poynting flux on the surface 2. The Poynting flux can be 
computed using the fact that £ is actually equal to the surface- 
tangential fields: £ = (n x H; — nxE) where n is the out- 
ward unit-normal vector. It follows that the integrated flux 
is-|Re^(ExH)-n= \ Re(£ 2 , 0°) (equivalent to the 
power exerted on the surface currents by the total field, with 
an additional 1/2 factor from a subtlety of evaluating the fields 
exactly on the surface [20]). Hence, 

$ s = I Re (^,/) 2 = i R e(C 2 ,0 2 ) 2 = -iRee-rV 2 ) 

where we used the continuity of and (j) 2 . Substituting £ 2 = 
J2 n x n@n an d recalling the definition (3) of G 2 , we obtain 
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via straightforward algebraic manipulations. 

Now, to obtain = (<E> S ) we must ensemble-average (• • • ) 
over all sources a 1 , and this corresponds to computing the ma- 
trix C — (ss*), which is only nonzero in its upper-left block 
C 1 = (s 1 s 1 *). Such a Hermitian matrix is completely de- 
termined by the values of x 1 *S(C 1 ) T Sx 1 for all vectors x 1 , 
where we have inserted the sign-flip matrices S and the trans- 
position for later convenience, and by study of this expression 
we will find that C 1 has a simple physical meaning. To begin 
with, we write £ x = to obtain: 
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where we have integrated over all possible dipole positions. 
The current-current correlation function (<j 1 (y)<j 1 (y / ) T ) = 
^S(y — y')w Im x is given by the fluctuation-dissipation the- 
orem [18], where we have factored out a 0(cj, T 1 ) term into 
Eq. (1) and where Imx denotes the imaginary part of the 
6x6 material susceptibility (whose diagonal blocks are Im e 
and Im //), related to material absorption (or the conductivity 
oo Im x). This eliminates one of the integrals, leaving 

£ J e (x')*Sn (x', y) [co Im x (y)] T 1 (x, yfS^ 1 (x). 
If we now employ reciprocity (from above), we can write 

where cf) 1 = T 1 * t; 1 is the field due to the surface current 
where the commuted S can be used to simplify the remaining 
term ^ 1 (x)*5T 1 (x, y)S = [r x (x, y)^(x)]*, assuming that 
S commutes with Im x (true unless there is a bi-anisotropic 
susceptibility, which breaks reciprocity). Finally, we obtain: 
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But ^ 1 *(cjlmx)0 1 = jj> Re f^ 1 * (-iuox^ 1 )] is exactly the 
time-average power density dissipated in the interior of ob- 
ject 1 by the field (j) 1 produced by £ x , since —iuox^ 1 is a 
bound-current density. 

Computing the interior dissipated power from an arbi- 
trary surface current is somewhat complicated, but mat- 
ters here simplify considerably because the C matrix is 
never used by itself — it is only used in the trace expres- 
sion $ = -|Tr[CW*(symG 2 )iy] = -\Tt[- - -} T = 
-\ Ty[SC t SW(symG 2 )W% by reciprocity as in Eq. (5). 
From the Cholesky factorization symG 2 = — U 2 *U 2 , this 
becomes \ Tr[X* SC T SX\, where X = WU 2 * are the "cur- 
rents" due to "sources" represented by the columns of £7 2 *, 
which are all of the form [0; s 2 ] (corresponding to sources 
in object 2 only). So, effectively, S (C 1 ) 7 ^ S is only used 
to evaluate the power dissipated in object 1 from sources in 
object 2, and by the same Poynting-theorem reasoning from 
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above, it follows that S (G 1 ) 
G 1 ~ ^>ymS{G 1 ) T S = - 1 symG 1 by the symmetry 
of G 1 , and Eq. (2) follows. 

It is also interesting to consider the spatial distribution of 
the Poynting-flux pattern, which can be obtained easily be- 
cause, as explained above, \ Re[£ 2 (x)*(/> 2 (x)] is exactly the 
inward Poynting flux at a point x on surface 2. It follows that 
the mean contribution <£ 2 of a basis function (5 T n to is 
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FIG. 1 : Flux spectra $ of isolated (d — > oo) and interacting (d = R) 
gold cylinders/spheres (solid/hollow circles) of length L and radii 
R = 0.2/xm, normalized by their corresponding surface areas A, 
computed via Eq. 2. Solid lines show <E> computed via the semi- 
analytical formulas of [10, 21]. 



where e 2 is the unit vector corresponding to the (3 2 compo- 
nent. This further simplifies to <£ 2 = F 2 n , where 
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Note that <I> = TrF 2 . Similarly, by swapping 1 <-> 2 we 
obtain a matrix F 1 such that (fr 1 = F^ n is the contribution 
of /3 1 to the flux on surface 1. In the case of BEM with the 
standard RWG basis [17], fi r n is localized around one edge of 
a triangular surface mesh, so the flux contribution of a single 
triangular panel can be computed from the sum of F^ n /2 from 
the edges of that triangle. 

For a single object 1 in medium 0, the emissivity of the ob- 
ject is the flux <£° of random sources in 1 into [6]. Follow- 
ing the derivation above, the flux into is —\ Re = 
~ \ (C 1 ? T° ★ The rest of the derivation is essentially un- 
changed except that W = (G 1 + G ) -1 since there is no 
second surface. Hence, we obtain 

$° = ^- Tr [(symG 1 ) W* (symG ) W] , (8) 

which again is invariant under 1 <H> interchange from the 
reciprocity relations (Kirchhoff 's law). 

Results: Figure 1 shows the flux spectrum $ for vari- 
ous configurations of gold spheres and cylinders (of radii 
R = 0.2/im and varying lengths L), as a function of fre- 
quency R/X. (<£ is normalized by the surface area A of each 
object to make comparisons easier. At these wavelengths, R 
is several times the skin depth S = c/^feuj, which means that 
most of the radiation is coming from sources near the sur- 
face [21].) Our results for isolated and interacting spheres 
(red hollow circles) agree with previous results based on semi- 
analytical formulas [10, 21] (solid lines). In addition, Fig. 1 
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FIG. 2: Flux spectra $ of isolated and interacting/interlocked 
spheres/rings (solid/hollow circles) of radii R — 1/im, normal- 
ized by their corresponding surface areas A, computed via Eq. 2. 
Solid lines denote $ as computed via the semi-analytical formulas 
of [10, 21]. Insets show the spatial distribution of surface flux pat- 
tern at particular frequencies (right colorbar). 



spection of the spatial distribution of flux pattern on the sur- 
face of the objects, which we compute via Eq. 7 and show as 
insets in Fig. 2, for both rings and spheres. As expected, at 
large wavelengths X^> R, near-field effects dominate and the 
flux pattern peaks in regions of nearest surfaces. However, for 
A ~ R, the sphere-sphere pattern does not change qualita- 
tively while the ring-ring pattern exhibits resonance patterns 
characterized by nodes and peaks distributed along the ring. 
(Interestingly, the flux pattern of the first resonance is peaked 
away from the nearest surfaces.) Away from these resonances, 
the ring emissivity is smaller: for A <C R (not shown), is 
well described by the Stephan-Boltzmann law, and the ratio 
of their emissivities is given by the ratio of their surface areas 
« 0.3. A similar reduction occurs for A ^> R due to the ring's 
smaller polarizability. 

This work was supported by DARPA Contract No. 
N66001-09-1-2070-DOD and by the AFOSR Multidisci- 
plinary Research Program of the University Research Initia- 
tive (MURI) for Complex and Robust On-chip Nanophoton- 
ics, Grant No. FA9550-09- 1-0704. 



shows $ for isolated and interacting cylinders (solid circles) 
of various aspect ratios L/R; previous results based on semi- 
analytical methods (solid lines) were limited to the infinite 
case L/R —> oo [21]. For L/R « 1 (not shown), correspond- 
ing to nearly-isotropic cylinders, is only slightly larger than 
that of an isolated sphere due to the small but non-negligible 
volume contribution to <£. As L/R increases, increases 
over all A, and converges towards the L — >> oo limit (black 
solid line) as A — >> 0, albeit slowly. Interestingly, ^> $oo 
at particular wavelengths, a consequence of geometrical reso- 
nances that are absent in the infinite case. (Away from these 
resonances, $ clearly straddles the L — » oo result so long as 
A < L.) For interacting cylinders, in addition to the expected 
near-field enhancement at large A, one also finds significant 
resonant peaks at A < L. 

Equation 2 can be exploited to obtain <1> in an even more 
complicated geometry, where the topology makes it difficult 
to distinguish the incoming and outgoing waves of other for- 
mulations [11-13]. Figure 2 shows <1> for isolated and in- 
terlocked gold rings (solid circles), of inner and outer radii 
r = 0.7/im and R = 1/im, respectively, and thickness 
h = 0.1/im. For comparison, we also show the correspond- 
ing for isolated and interacting spheres of radii R (open 
circles). As in the case of finite cylinders, the rings exhibit 
orders of magnitude enhancement in at particular A, corre- 
sponding to azimuthal resonances — the first of which is the 
m = mode at A « 2ttR. Interestingly, despite its smaller 
surface area and volume, the absolute (unnormalized) <E> of 
the isolated ring is « 4.5 times larger than that of the sphere 
at the fundamental resonance. The geometrical origin of this 
resonance enhancement becomes even more apparent upon in- 
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